      PROGRAM TRACKING

      common/bl1/ dx,dy
      common/bl2/ vector(2,512,512),velin(2),points(2,64000)

      character*16 ARCH1,ARCH2,ARCH3
      character*16 TIEMPO
      character TITULO*80
      real rmv(64000)

CCC----------------------------------------------------------
CCC  	Lectura de los Datos
CCC----------------------------------------------------------
      open(10,FILE='DatosPT',STATUS='UNKNOWN')
      read(10,*)
C      read(10,*) icont
C      read(10,*)

C
C Longitudes de la cavidad
C
      read(10,*) flx, fly
      read(10,*)
C
C Numero de nodos
C
      read(10,*) NCI, NCJ
      read(10,*)
C
C Numero de trayectorias a calcular
C
      read(10,*) Npoints
      read(10,*)
C
C Numero de archivo de velocidades inicial 
C
      read(10,*) init
      read(10,*)
C
C Numero de archivo de velicidades final
C
      read(10,*) NMAX
      read(10,*)
C
C Salto en la lectura de los archivos
C
      read(10,*) NPR
      read(10,*)
C
C Incremento de tiempo con el que se calcularon los campos de velocidad
C
      read(10,*) DT
      read(10,*)
C
C Numero de ciclos a calcular
C
      read(10,*) NREP
      read(10,*)
C
C Frecuencia de impresion
C
      read(10,*) NFREQ
      close(10)

      write(*,*) " "
      write(*,*) " "
      write(*,*) " ---------------------"
      write(*,*) " | PARTICLE TRACKING |"
      write(*,*) " ---------------------"
      write(*,*) "  "
      write(*,*) "  "
      write(*,*) " Numero de trayectorias = ", Npoints
      write(*,*) " Inicio = ", init
      write(*,*) " Final = ", NMAX
      write(*,*) " Salto entre archivos = ", NPR
      write(*,*) " DT = ", DT
      write(*,*) " NCI = ", NCI
      write(*,*) " NCJ = ", NCJ
      write(*,*) " Longitud x = ", flx
      write(*,*) " Longitud y = ", fly
      write(*,*) " Numero de ciclos = ", NREP
      write(*,*) " Frecuencia de impresion = ", NFREQ

CCC----------------------------------------------------------
CCC Incrementos de la malla de la cavidad
CCC----------------------------------------------------------
      dx = flx/float(NCI-1)
      dy = fly/float(NCJ-1)

CCC----------------------------------------------------------
CCC Posiciones de las particulas
CCC----------------------------------------------------------

      if (NCONTINU .eq. 1) then
CCC----------------------------------------------------------
CCC Lee los puntos de un archivo
CCC----------------------------------------------------------
         IU = MOD(init,10)
         ITOT = (ncont - IU)/10
         ID = MOD(ITOT,10)
         ITOT = (ITOT - ID)/10
         IC = MOD(ITOT,10)
         ITOT = (ITOT - IC)/10
         IM = MOD(ITOT,10)
         ITOT = (ITOT - IM)/10
         IX = MOD(ITOT,10)
         WRITE(TIEMPO(1:1),3000) IX
         WRITE(TIEMPO(2:2),3000) IM
         WRITE(TIEMPO(3:3),3000) IC
         WRITE(TIEMPO(4:4),3000) ID
         WRITE(TIEMPO(5:5),3000) IU
         ARCH1(6:10) = TIEMPO
         open(50,FILE=ARCH1,STATUS='UNKNOWN')    
         do i = 1,Npoints
            read(50,3040) points(1,i), points(2,i),velocity
         enddo
         nt = init
         close(50)
CCC----------------------------------------------------------
CCC Construye la linea de puntos
CCC----------------------------------------------------------      
      else
         dxpoint = (flx - 2 * dx) / (Npoints - 1)
         write(*,*) " dxpoint = ", dxpoint
         write(*,*) " dx = ", dx
         open(50,FILE="line.00000",STATUS='UNKNOWN')          
         do i = 1,Npoints
            points(1,i) = dx + dxpoint * (i-1)
            points(2,i) = fly * 0.5
            write(50,3040) points(1,i),points(2,i),0.0 
         enddo
         close(50)
         nt = 0
      endif
      
CCC----------------------------------------------------------
CCC Prefijo de los archivos de velocidad y de las geometrias
CCC----------------------------------------------------------
      ARCH2 = 'velc.'
      ARCH3 = 'line.' 

      nsteps = NMAX
C
C incremento de tiempo que se usara en el calculo
C i.e. en el metodo de Runge-Kutta
C
      DT1 = DT*NPR

CCC********************** Ciclo principal ***********************
CCC----------------------------------------------------------
CCC El primer loop es para repetir el proceso nrep veces
CCC----------------------------------------------------------
      do irep = 1,nrep
         do ncont=init,nsteps,npr
CCC----------------------------------------------------------
CCC Genera el nombre del archivo que se va a leer (velc.)
CCC----------------------------------------------------------
            IU = MOD(ncont,10)
            ITOT = (ncont - IU)/10
            ID = MOD(ITOT,10)
            ITOT = (ITOT - ID)/10
            IC = MOD(ITOT,10)
            ITOT = (ITOT - IC)/10
            IM = MOD(ITOT,10)
            ITOT = (ITOT - IM)/10
            IX = MOD(ITOT,10)
            WRITE(TIEMPO(1:1),3000) IX
            WRITE(TIEMPO(2:2),3000) IM
            WRITE(TIEMPO(3:3),3000) IC
            WRITE(TIEMPO(4:4),3000) ID
            WRITE(TIEMPO(5:5),3000) IU
            ARCH2(6:10) = TIEMPO
            open(30,FILE=ARCH2,STATUS='UNKNOWN')
            read(30,*) (((vector(IJ,I,J),IJ=1,2),I=1,NCI),J=1,NCJ)

CCC----------------------------------------------------------
CCC     RUNGE-KUTTA
CCC----------------------------------------------------------
            do iw=1,NPoints
               xs = points(1,iw)
               ys = points(2,iw)

               call INTER(iw)
               a1 = velin(1)*DT1
               a2 = velin(2)*DT1
               points(1,iw) = xs + a1/2.0
               points(2,iw) = ys + a2/2.0              

               call INTER(iw)
               b1 = velin(1)*DT1
               b2 = velin(2)*DT1
               points(1,iw) = xs + b1/2.0
               points(2,iw) = ys + b2/2.0
               
               call INTER(iw)
               c1 = velin(1)*DT1
               c2 = velin(2)*DT1
               points(1,iw) = xs + c1
               points(2,iw) = ys + c2

               call INTER(iw)
               d1 = velin(1)*DT1
               d2 = velin(2)*DT1
               
               xs = xs + (a1 + 2*b1 + 2*c1 + d1)/6.0
               ys = ys + (a2 + 2*b2 + 2*c2 + d2)/6.0
               if( xs .GE. 1) xs = flx
               if( ys .GE. 1) ys = fly
               if( xs .LE. 0) xs = 0.0
               if( ys .LE. 0) ys = 0.0
               points(1,iw) = xs
               points(2,iw) = ys

               rmv(iw) = SQRT(velin(1)*velin(1) + velin(2)*velin(2))
            enddo

CCC----------------------------------------------------------
CCC Genera el nombre del archivo que se va a guardar (line.)
CCC----------------------------------------------------------
            nt = nt + 1
            if( MOD(nt,NFREQ) .eq. 0 ) then
               IU = MOD(nt,10)
               ITOT = (nt - IU)/10
               ID = MOD(ITOT,10)
               ITOT = (ITOT - ID)/10
               IC = MOD(ITOT,10)
               ITOT = (ITOT - IC)/10
               IM = MOD(ITOT,10)
               ITOT = (ITOT - IM)/10
               IX = MOD(ITOT,10)
               WRITE(TIEMPO(1:1),3000) IX
               WRITE(TIEMPO(2:2),3000) IM
               WRITE(TIEMPO(3:3),3000) IC
               WRITE(TIEMPO(4:4),3000) ID
               WRITE(TIEMPO(5:5),3000) IU
               ARCH3(6:10) = TIEMPO
               write(*,*) "Storing file : ", ARCH3
               open(40,FILE=ARCH3,STATUS='UNKNOWN')            
               
               do iw=1,NPoints
ccc                  write(40,3040) points(1,iw),points(2,iw),rmv(iw) 
                  write(40,3040) points(1,iw),points(2,iw) 
               enddo
               close(40)
            endif

            write(*,3015) TIEMPO
            
            close(30)
         enddo        

      enddo      
CCC********************************************************
	
 3000 format(i1)
 3010 format(a80)
 3015 format(a5)
 3030 format(6e12.5)
ccc 3040 format(3e12.5)
 3040 format(2e12.5)
      
      write(*,*) "Todo bien !!!"

      stop
      end

      subroutine INTER(ip)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       Subrutina para la interpolacion usando una serie
C       de Taylor a primer orden. Recibe como parametro la
C       posicion de la particula en el instante actual y
C       calcula la velocidad en esa posicion.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                 
      common/bl1/ dx,dy
      common/bl2/ vector(2,512,512),velin(2),points(2,64000)

        i = NINT(points(1,ip)/dx)
        j = NINT(points(2,ip)/dy)
        a = points(1,ip)/dx
        b = points(2,ip)/dy
        
        alpha = a - i 
        beta = b - j

        velin(1) = ( vector(1,i,j) * (1.0 - alpha) +
     &               vector(1,i+1,j) * alpha) * (1.0 - beta) +
     &             ( vector(1,i,j+1) * (1.0 - alpha) +
     &               vector(1,i+1,j+1) * alpha ) * beta

        velin(2) = ( vector(2,i,j) * (1.0 - alpha) +
     &               vector(2,i+1,j) * alpha) * (1.0 - beta) +
     &             ( vector(2,i,j+1) * (1.0 - alpha) +
     &               vector(2,i+1,j+1) * alpha ) * beta

        return
        end
